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The  use  of  the  Global  Positioning  System  (GPS)  as  a 
navigation  aid  for  aircraft  attempting  to  locate  a  ground-based 
electromagnetic  energy  emitter  is  studied.  In  particular,  the 
satellite  geometry  which  yields  minimum  errors  in  the  emitter 
location  estimation  for  four  different  satellite  availability 
cases  is  explored.  This  geometry,  in  general,  is  not  the  same 
as  that  which  yields  minimum  aircraft  navigation  errors. 

Satellite  selection  criteria  are  identified  and  serve  as  a  basis 
for  selection  algorithm  development. 

The  research  shows  that  emitter  location  errors  can  be 
significantly  reduced  by  selecting  satellites  based  on  the 
criteria  presented  in  this  study.  Three  satellite  performance 
is  found  to  be  nearly  as  good  as  that  obtained  using  four 
satellites  and,  for  the  two  satellite  case,  emitter  location  is 
still  better  for  some  period  of  time  than  that  obtained  using 
four  satellites  selected  to  minimize  Geometric  Dilution  of 
Precision  (GDOP) . 


OPTIMAL  SELECTION  OF  GLOBAL  POSITIONING  SYSTEM  SET 
TO  MINIMIZE  EMITTER  LOCATION  ERRORS 

I.  INTRODUCTION 

1.1  BACKGROUND 

The  use  of  aircraft  equipped  with  electro-magnetic  energy 
receivers  to  determine  the  location  of  an  emitter  has  many 
military,  as  well  as  civilian,  applications.  Missions  such  as 
strategic  and  tactical  reconnaissance,  as  well  as  airborne 
search  and  rescue,  require  accurate  determination  of  an  emitter 
location.  The  difficulty  faced  by  these  airborne  collectors  is 
that  in  order  to  determine  an  emitter  location  precisely,  they 
must  establish  their  own  position  very  accurately.  This 
difficulty  is  compounded  by  the  fact  that,  for  many 
applications,  a  global  positioning  capability  is  required. 
Although  many  navigation  systems  such  as  LORAN  and  OMEGA  provide 
substantial  capability,  the  Global  Positioning  System  (GPS) 
provides  continuous  global  coverage  with  much  greater  accuracy. 
However,  in  order  to  obtain  highly  accurate  observer  position 
information,  it  is  necessary  that  the  airborne  platform  select 
the  appropriate  satellites  from  the  available  constellation  of 
satellites  which  will  yield  optimum  emitter  location  solution 
geometry.  Ideally,  the  satellite  selection  criterion  should 
optimize  the  accuracy  of  the  emitter  location  estimation. 


1.2  PROBLEM  AND  SCOPE 

The  problem  is  to  determine  the  positions  of  three  airborne 
collectors  and  four  GPS  satellites  such  that  errors  in  the 
estimate  of  emitter  location  will  be  minimized.  In  this  case, 
optimum  position  is  defined  as  that  position  which  results  in  a 
minimum  value  for  the  emitter  location  estimate  mean  squared 
miss  distance  (MSMD) ,  or  alternatively,  the  circular  error 
probable  (CEP) .  The  MSMD  is  a  scalar  cost  function  of  the 
emitter  location  errors,  which  in  turn  is  a  function  of  the 
satellite  position  geometry  expressed  in  the  Kalman  filter 
measurement  equation.  The  local  minimum  CEP  cost  value  is 
obtained  using  a  steepest  descent  gradient  search  algorithm. 

This  study  examines  the  impact  of  both  satellite  positions 
and  collector  positions  on  emitter  location  errors.  It  also 
investigates  emitter  location  accuracy  after  degrading  to  three 
and  two  satellite  operation.  Additionally,  an  algorithm  is 
developed  which  selects  the  optimum  satellite  configuration  for 
minimum  emitter  location  errors.  Note  that  the  goal  is  not  to 
develop  the  "best  possible"  emitter  location  system,  but  to 
determine  the  optimal  satellite  geometry  which  minimizes  emitter 
location  error.  Therefore,  a  highly  elaborate  truth  model  is 
not  required. 


1.3  SUMMARY  OF  CURRENT  KNOWLEDGE 

The  NAVSTAR  Global  Positioning  System  is  a  space-based 
radio-positioning  navigation  system  which  provides  highly 
accurate  three-dimensional  position  and  velocity  information  on 
a  global  basis  to  a  very  large  number  of  users  (9:146).  In 
general,  at  least  four  satellites  are  required  to  solve  the  four 
time-difference-of-arrival  equations  for  the  three  dimensional 
position  and  the  user  receiver  clock  error  (1:85).  Since  the 
GPS  generally  allows  the  user  to  view  six  or  more  satellites  at 
any  given  time,  it  is  necessary  to  select  the  combination  of 
satellites  that  will  give  the  most  accurate  position  information 
(3:8).  The  effect  that  satellite  geometry  has  on  ranging 
accuracy  can  be  expressed  in  terms  of  Geometric  Dilution  of 
Precision  (GDOP)  (3:8).  An  algorithm  which  maximizes  navigation 
accuracy  selects  the  combination  of  four  satellites  that  yields 
the  smallest  value  of  GDOP  (3:8).  It  should  be  noted  that 
although  the  minimum  GDOP  criteria  results  in  minimum  user 
position  errors,  it  "...  does  not  generally  result  in  minimum 
emitter  location  errors"  (4:64).  This  is  due  to  two  factors. 

One  is  that  vertical  collector  position  errors  contribute 
negligibly  to  emitter  location  errors  since  the  emitter  is 
assumed  to  be  on  the  surface  of  the  earth.  The  second  is  that 
the  collector  position  errors  are  mapped  through  the  nonlinear 
hyperbolic  functions  to  the  emitter  location  errors  (10) .  The 


use  of  joint  estimation  of  emitter  position  with  observer 


position  produces  emitter  location  errors  which  are 
significantly  smaller  than  estimates  based  on  satellite 
selection  using  minimum  GDOP  criteria  (4:93).  It  is  also  noted 
that  the  joint  estimator  significantly  outperforms  the  two  stage 
estimation  process  where  the  observer  navigation  problem  is 
solved  first,  followed  by  solution  for  the  emitter  position. 

The  price  which  is  paid  for  the  gains  obtained  by  using  a  joint 
estimator  is  an  increase  in  computational  loading.  During  his 
research  at  the  Massachusetts  Institute  of  Technology, 
Lewantowicz  (4)  investigated  system  performance  when  using  only 
three  satellites  and  found  that  emitter  location  accuracy  was 
very  nearly  the  same  as  that  obtained  using  four  satellites. 
Additionally,  he  demonstrated  that  degraded  operation  using  only 
two  satellites  and  a  precise  clock  yields  an  acceptably  accurate 
emitter  location  solution. 

A  number  of  assumptions  necessary  to  define  and  adequately 
limit  the  problem  are  presented  in  the  next  section. 

1.4  ASSUMPTIONS 

The  following  assumptions  are  made  at  the  outset  of  the 
problem: 

1.  There  are  three  airborne  collectors  operating 
simultaneously . 

2.  Each  collector  platform  has  a  high  grade  Inertial  Navigation 
System  (INS)  with  barometric  altimeter  data  available.  The 


external  altitude  measurement  is  available  to  stabilize  the 
inherently  unstable  vertical  channel  of  the  INS. 

3.  Each  collector  has  measurements  available  from  four  GPS 
satellites.  Further,  each  observer  uses  the  same  four 
satellites  for  updates.  This  is  a  reasonable  assumption  since 
the  three  platforms  are  operating  in  the  same  geographic  area 
and  will  "see"  essentially  the  same  constellation  of  available 
satellites. 

4.  The  GPS  satellites  selected  should  lie  at  or  above  five 
degrees  elevation  from  the  local  horizon  measured  from  the 
center  of  the  collector  geometry.  Lower  elevation  angles 
increase  ranging  errors  due  to  propagation  effects. 

5.  The  collectors  use  a  passive  emitter  locating  method  known 
as  the  hyperbolic  location  system  (10) . 

6.  A  constrained  emitter-collector  geometry  is  assumed.  This 
assumption  narrows  the  problem  to  that  where  the  approximate 
emitter  location  in  general  is  initially  known  or  crudely 
measured. 

1.5  GENERAL  APPROACH 

Parameter  optimization  is  used  to  solve  a  set  of  emitter 
location  problems  assuming  various  satellite  availability 
constraints  as  outlined  by  the  following. 


! 


OPTIMUM  SATELLITE  POSITION 


During  this  portion  of  the  study  the  four  GPS  satellite 
positions  are  selected  from  among  an  infinite  set  on  the  orbital 
sphere  using  the  cost  gradient  search.  The  collectors  executed 
their  predetermined  orbits.  The  cost  is  formulated  as  a 
function  of  GPS  satellite  positions  and  is  computed  using  the 
Kalman  filter  covariance  matrix  function.  The  filter  covariance 
matrix  is  a  function  of  the  measurement  observation  matrix, 
which  in  turn  depends  on  the  satellite  positions.  Thus,  for  a 
steady-state  Kalman  filter  solution,  the  cost  function  depends 
only  on  the  satellite  positions.  Therefore,  the  gradient  vector 
of  this  scalar  cost  function  indicates  the  magnitude  and 
"direction  of  movement"  for  each  satellite  at  update  time.  The 
navigation  Kalman  filter  is  allowed  to  reach  steady  state  from 
the  same  initial  conditions  at  each  cost  computation  point 
before  satellite  "movements"  are  computed.  It  is  hypothesized 
that  information  retained  by  iterating  the  filter  only  one  time 
between  satellite  movements  without  reinitializing  the  Kalman 
filter  may  have  produced  overly-optimistic  CEPs  as  obtained  by 
Lewantowicz  in  his  research  (4) .  His  motivation  for  iterating 
only  one  time  is  the  computational  savings  in  the  minimum  cost 
search. 

THREE  SATELLITE  PERFORMANCE 

Of  the  four  satellites,  three  are  moved  to  optimum 
positions  for  this  portion  of  the  study,  while  the  fourth 
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satellite  remains  fixed  overhead  to  simulate  the  availability  of 
a  precise  clock  onboard  each  collector  platform.  A  second 
analysis  is  accomplished  using  information  from  only  three 
satellites  to  determine  the  validity  of  the  hypothesis  that  the 
overhead  satellite  can  be  used  to  simulate  a  precise  clock  in  a 
three  satellite  problem.  A  comparison  is  made  between  the 
results  obtained  using  these  two  approaches. 


TWO  SATELLITE  PERFORMANCE 


This  portion  of  the  study  is  very  similar  to  the  three 
satellite  case.  Given  Kalman  filter  estimates  of  user  clock 
bias  and  user  clock  drift  obtained  using  measurements  from  four 
satellites,  the  number  of  visible  satellites  is  reduced  to  only 
two  and  the  satellite  positions  are  again  optimized. 

Comparisons  are  made  with  previous  simulations  to  determine  the 
feasibility  of  two  satellite  operation. 


DEVELOPMENT  OF  A  SELECTION  ALGORITHM 

The  results  of  the  optimum  satellite  position  analysis  are 
then  used  to  construct  a  selection  algorithm.  The  behavior  of 
the  scalar  cost  as  a  function  of  satellite  geometries  are 
studied  to  determine  the  characteristics  of  an  "optimum 
satellite  geometry"  and  selection  criteria  are  established. 


1.6  THESIS 


The  following  chapters  describe  the  problem  structure,  the 
mechanics  used  in  solving  the  problem,  and  the  results  obtained. 

Chapter  2  develops  the  error  dynamics  and  measurement 
models  used  by  the  emitter  locating  system.  The  observer 
navigation  error  state  vector  is  augmented  with  the  error  states 
associated  with  the  emitter  measurement  process  to  obtain  the 
system  error  state  vector.  Measurement  models  are  formulated 
for  observer  navigation  measurements,  as  well  as  emitter 
measurements,  and  are  linearized  for  use  in  linear  estimation. 

Chapter  3  builds  the  structure  of  the  discrete-time 
extended  Kalman  filter  which  is  used  as  the  estimator  for  this 
study.  The  cost  function,  based  upon  emitter  location  mean 
squared  miss  distance,  is  developed  as  a  function  of  the  GPS 
satellite  positions.  The  emitter  location  CEP  is  also  computed. 
Finally,  the  minimum  cost  search  technique  used  in  this  study, 
the  steepest  descent  weighted  gradient  algorithm,  is  described. 

Chapter  4  presents  the  results  obtained  when  the  GPS 
satellites  are  moved  to  optimum  positions,  minimizing  the  error 
in  the  emitter  location  estimate.  The  three  satellite 
performance  is  shown  to  be  essentially  the  same  as  that  obtained 
using  four  satellites.  Further,  the  use  of  only  two  satellites 


and  filter  estimates  for  user  clock  parameters  yields  reasonable 
emitter  location  performance.  Satellite  selection  criteria  are 
established  based  on  the  results  obtained  in  the  four  satellite 
optimization  cases. 

Chapter  5  presents  conclusions  and  recommendations  arising 


from  this  study. 


II.  ERROR  DYNAMICS  AND  MEASUREMENT  MODELS 

2.1  INTRODUCTION 


In  general,  locating  an  electromagnetic  emitter  requires  a 
method  for  determining  the  observer's  position,  as  well  as  a 
method  for  distinguishing  signals  intercepted  from  an  emitter. 
For  the  simplified  case  of  a  stationary  observer,  the  observer 
position  can  be  very  accurately  determined  by  means  of  a  precise 
survey.  However,  for  such  applications  as  airborne  search  and 
rescue  operations,  the  observers  must  move  in  order  to  provide 
the  required  coverage  of  the  vast  areas  involved.  The  problem 
of  determining  the  position  of  these  moving  observers  requires  a 
navigation  system  capable  of  providing  highly  accurate  position 
information  on  a  global  basis. 

In  this  study,  the  system  used  by  each  observer  to  provide 
this  high  quality  position  data  is  a  high  grade  inertial 
navigation  system  (INS)  which  is  updated  by  measurements  from 
four  GPS  space  vehicles  (SV)  and  stabilized  by  measurements 
from  a  barometric  altimeter.  As  pointed  out  in  Section  1.4,  the 
GPS  satellites  selected  lie  at  or  above  five  degrees  elevation 
from  the  local  horizon  and  each  observer  uses  the  same  four 
satellites  for  updates.  Note  that  the  following  models  and 
derivations  correspond  to  those  used  and  developed  by 
Lewantowicz  in  his  thesis  (4) . 


The  method  used  to  analyze  the  performance  of  the  emitter 
location  system  is  a  linearized  error  covariance  analysis  of  the 
estimated  errors.  To  perform  this  analysis,  the  state  vector 
representing  the  actual  system  state  is  given  by  j<a(t)  and  the 
estimate  of  the  state  vector  at  the  same  time  t  is  given  by 

A 

2S(t)  .  The  error  state  vector,  x(t),  can  then  be  defined  as 

2£(t)  -  xa(t)  -  x(t)  (2-1) 

which  satisfies  in  general  the  vector  differential  equation 

X(t)  -  l(x,t)  +  G  w(t)  (2-2) 

X(0)  -  x0 

where  f  is  a  time-varying  vector  valued  function  of  x(t) 
describing  the  error  state  dynamics.  The  matrix  G  is  the  time 
invariant  noise  distribution  matrix,  and  w  is  an  independent 
zero-mean  white  Gaussian  noise  process  with  covariance  kernel 

E  [  W  ( t )  WT  { t  +  t)]  =  Q(t)  Mr)  (2-3) 

where  E  is  the  expectation  operator,  Q  is  a  diagonal  matrix,  anc 
is  the  Dirac  delta  function. 

The  measurement  process  for  the  estimator  is  well  modelled 


where  &  is  a  time-varying  vector  valued  function  of  the  error 
state  x.  The  measurement  noise  v(t)  is  an  independent  zero-mean 
white  Gaussian  noise  process  with  covariance  kernel 


E[v(t)yT(t  +  t ) ]  *  R  5 (  t )  (2-5) 

where  R  is  a  diagonal  matrix. 

Given  this  system  model,  an  extended  Kalman  filter  is  used 
to  estimate  navigation  errors.  This  filter  is  discussed  in 
Section  3.2.  System  dynamics  and  measurement  models  for 
navigation  and  emitter  location  are  discussed  in  Sections  2.2 
through  2.4. 

2.2  OBSERVER  NAVIGATION  SYSTEM 

A  number  of  systems  are  available  which  can  provide 
navigation  data  to  airborne  observers.  These  include  time- 
dif ference-of-arrival  (TDOA)  type  systems  such  as  LORAN  and 
OMEGA  and  range  measurement  type  systems  using  Distance 
Measuring  Equipment  (DME)  measurements  to  surveyed  stations  on 
the  ground.  These  systems  are,  in  general,  limited  by 
combinations  of  accuracy,  availability  or  areas  of  coverage.  A 
performance  improvement  over  these  systems  is  provided  by  GPS. 


The  NAVSTAR  GPS,  when  fully  operational,  will  consist  of  a 
constellation  of  18  SVs,  3  in  each  of  6  orbital  planes  (9:146). 
The  GPS  has  potential  for  providing  highly  accurate  three- 
dimensional  position  and  velocity  information  along  with 
Coordinated  Universal  Time  (UTC)  to  a  very  large  number  of 
suitably  equipped  users.  The  orientation  of  the  satellite 
orbits  generally  allows  the  observer  to  view  six  or  more 
satellites  at  any  given  time,  thus  providing  adequate  navigation 
information  on  a  global  basis  (3:18).  The  GPS  SV  transmits  an 
encoded  navigation  message  from  which  the  receiver  can  determine 
the  pseudo-range  to  the  SV  with  a  very  high  degree  of  accuracy. 
The  locus  of  points  for  each  pseudo-range  measurement  describes 
a  sphere  whose  center  is  at  the  SV.  The  point  where  three 
spheres  intersect  provides  the  observer  with  position 
information  in  three  dimensions.  Since  the  observer  generally 
uses  a  fairly  inaccurate  crystal  clock,  a  fourth  pseudo-range 
measurement  is  used  to  determine  the  user  clock  bias.  Most  of 
the  error  sources  associated  with  the  GPS  can  be  effectively 
minimized  by  error  modelling  and  the  use  of  Kalman  filtering  to 
estimate  these  errors.  As  a  result,  the  GPS  is  capable  of 
providing  much  more  accurate  position  information  than  most 
other  navigation  systems.  Because  of  this  high  degree  of 
accuracy  and  the  global  positioning  capability  of  the  GPS,  a 
high-grade  INS  is  combined  with  the  GPS  to  form  the  navigator 
which  is  used  for  the  emitter  locating  system  in  this  study. 


2.2.1  THE  INERTIAL  NAVIGATION  SYSTEM 


The  INS  system  errors  are  modelled  using  16  states,  forming 
the  state  vectors  x1(  x2,  and  x3.  Note  that  there  are  three  INS 
error  state  vectors  since  each  of  the  three  airborne  observers 
has  its  own  INS.  Each  of  these  state  vectors  is  defined  as 


follows  (4) : 


*1 


East  position  error 
North  position  error 
Vertical  position  error 
East  velocity  error 
North  velocity  error 
Vertical  velocity  error 
Roll  error  (east  axis) 

Roll  error  (north  axis) 

Roll  error  (vertical  axis) 

Barometric  altimeter  error 

Barometric  sea  level  pressure  variation 

East  accelerometer  bias 

North  accelerometer  bias 

East  gyro  bias 

North  gyro  bias 

Vertical  gyro  bias 


(2-6) 


State  vectors  x2  and  £3  are  defined  in  the  same  manner.  It 


should  be  noted  that  the  roll,  pitch  and  yaw  data  in  each  body 


frame  are  transformed  to  the  East-North-Vertical  (ENV) 


navigation  frame  for  each  observer  (4:17). 


2.2.2  THE  GLOBAL  POSITIONING  SYSTEM 


The  measurements  available  to  update  the  inertial 
navigation  system  are  pseudo-range  measurements  to  the  GPS  SV 


where  R^j  is  the  pseudo-range  from  the  i-th  observer  to  the  j-th 

SV,  c  is  the  speed  of  light,  and  tij  is  the  travel  time  of  the 

signal.  Because  there  are  errors  inherent  in  the  system,  this 

pseudo-range  measurement  does  not  only  represent  the  actual  line 

of  sight  (LOS)  range  to  the  SV,  but  includes  several  errors. 

The  error  sources  modelled  for  this  study  are  uncalibrated 

propagation  errors  in  each  GPS  receiver  channel,  uncalibrated 

offset  and  drift  errors  in  each  observer  clock,  code  loop 

interchannel  biases,  and  LOS  errors.  These  errors  form  the 

state  vectors  x4  through  xg  described  by 

x4  12  states  for  receiver  propagation  errors 
X5  3  states  for  clock  offsets 

x6  3  states  for  clock  drifts 

x7  12  states  for  code  loop  interchannel  bias 
x8  4  states  for  LOS  biases  (l  per  SV) 

The  augmented  error  state  vector  is  formed  by  combining  all 
the  error  states  as 


TT'TTTTTT 

£-1'  X2,  X3,  X4;  x5,  X6,  X7,  X8] 


(2-8) 


The  resulting  error  state  column  vector  is  of  dimension  82. 

Although  models  of  higher  dimensionality  exist,  the  model 
described  is  adequate  for  this  study.  Such  assumptions  as 
straight,  level,  unaccelerated  flight  at  high  altitude  reduce 
the  effect  of  other  error  sources  allowing  lower  state  vector 


dimensionality. 


2.3  THE  EMITTER  MEASUREMENT  SYSTEM 

States  which  describe  the  errors  associated  with  the 

emitter  measurement  process  are  added  to  the  model .  These 

states  are  the  emitter  position,  emitter-signal  receiver  time 

delay  calibration  error  for  each  observer,  and  the  error  in 

modelling  tropospheric  delay  along  the  line  of  sight  from  each 

observer  to  the  emitter.  These  errors  form  the  state  vectors  x9 

through  x1:L  and  are  described  as 

x9  3  states  for  emitter  position  errors 
x10  3  states  for  receiver  calibration  errors 

x^  3  states  for  tropospheric  delay  errors 

The  total  augmented  error  state  vector,  x(t)  is  formed  by 
augmenting  state  vectors  x9  through  x1;L  to  Eqn.  (2-8) 

X  =  [XX,  X2,  X3,  X4,  X5,  X6,  X7,  Xg,  Xg,  xi0,  X^]  (2-9) 
This  final  error  state  vector  is  of  dimension  91. 

Note  that  this  formulation  of  the  error  state  vector  allows 
joint  estimation  of  emitter  position  with  observer  position  so 
that  all  available  information  is  processed  jointly  at  each 
computation  step.  The  result  of  this  is  that  the  solution  is 
optimal  in  the  minimum  mean-squared  error  (MMSE)  sense  (6) . 

Each  measurement  updates  both  the  navigation  position  and  the 
emitter  location.  As  a  result,  the  emitter  measurement  can  be 
used  to  reduce  navigation  errors.  The  consequence  of  choosing 


the  joint  estimation  scheme 
increase  in  required  computations.  For  this  study,  the 
increased  performance  of  the  joint  estimator  justified  the 
increased  computational  burden. 

With  the  models  for  the  error  state  dynamics  and  noise 
inputs  now  developed,  the  next  step  is  to  model  the  measurements 
for  the  navigation  update  and  the  emitter  update. 


2.4  OBSERVER  NAVIGATION  AND  EMITTER  MEASUREMENTS 


2.4.1  OBSERVER  NAVIGATION  MEASUREMENTS 

The  actual  navigation  measurement  is  the  pseudo-range  from 
the  observer  to  the  GPS  SV  defined  as  the  transit  time  of  the 
signal  scaled  by  the  speed  of  light.  Each  GPS  signal  carries  an 
encoded  navigation  message  containing  SV  ephemeris  data  which 
allows  the  position  of  the  SV  to  be  determined  very  accurately. 
With  this  information,  the  observer  can  solve  four  equations  in 
four  unknowns.  These  unknowns  are  the  three  position  components 
of  the  navigation  error  state  vector  and  the  user  clock  bias 
error. 

In  addition  to  the  inherent  system  errors  in  the  GPS,  other 


error  sources  affect  the  accuracy  of  the  pseudo-range 
measurement.  The  significant  errors  are  modelled  and  discussed 


in  this  section.  The  errors  induced  by  atmospheric  delay  result 
primarily  from  ionospheric  refraction.  The  change  in  apparent 
path  length  brought  about  by  ionospheric  refraction  can  be 
substantially  reduced  by  employing  dual  frequency  compensation 
since  each  SV  transmits  two  frequencies,  f1  and  f2.  The  two 
frequency  compensation  is  based  on  the  fact  that  the  ionospheric 
delay  8t  is  frequency  dependent 


8  t  = 


(2-10) 


where 

K  =  environmental  constant 
f  =  carrier  frequency 

Let  At1  denote  transit  time  measurement  at  frequency  f-^  and  write 
it  as  (4:23) 


Atj  =  At  +  +  r^ 


(2-11) 


where 


At  *  the  transit  time  of  the  signal  (uncorrupted) 
»s  ionospheric  delay  at  frequency  fi 
r-L  =  time  measurement  error  due  to  other  sources 


Forming  the  ratio 


Rewriting  Eqn.  (2-11)  for  frequency  f2 


At2  =  At  +  5t2  +  r2  (2-14) 

Next,  define  the  measurement  difference  between  frequency  f-^  and 
f2,  5  At  as 


SAt  =  At2  -  At2  (2-15) 

Solving  for  SAt  by  substituting  Eqn.  (2-11)  and  (2-14)  into 
(2-15)  yields 

SAt  =  At  +  5r1  +  rA  -  At  -  5t2  -  r2 

e  (2-16) 
=  8ti  -  8t2  +  rl  -  r2 

Substituting  (2-13)  for  St2  produces 


SAt  =  5t.  - 


T2  5t1  +  rl 


-  r- 


=  (1  -  Cl2)  5 T1 


+  r->  -  r- 


(2-17) 
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where 


o  = 


f2 


Writing  the  compensated  time  delay  measurement  at  frequency  f-^ 


'■S*' 

Atc  =  At,  - 


SAt 


-  n  2 


(2-11 


l-o 


Substituting  (2-17)  and  (2-11)  into  (2-18)  gives 


Atc  «  At  +  St 1  +  r1  - 


5  At 


l-o 


At  +  +  rjL  -  <$1^  - 


rl  -  r2 


l-o 


=  At  + 


ri  (l-o^)  -  r,  +  r- 


l-o 


-  At  - 


(2-1 


Note  that  the  ionospheric  delay  dependence  has  vanished,  but  r 
and  r->  are  still  undefined. 


.  iVi'j  I’t  I*. 


[">  t‘.  tV|V  !*»  I  .  I*.  iV.I’t 


vwr^wx¥nr"¥nmrTyT'| 


rl  -  yl  +  Atlc  +  V1 


r2  =  y2  +  At2c  +  v2 


(2-20) 


(2-21) 


where  )'1  and  Y2  are  the  GPS  code  loop  measurement  bias  errors  in 
the  two  receivers.  These  errors  are  modelled  as  elements  of 
state  vector  x4  while  Atc,  the  user  clock  bias  error,  is 
modelled  by  x5,  as  described  in  Section  2.2.2.  The  terms  vx  and 
v2  represent  zero  mean,  white  Gaussian  measurement  noise. 


Substituting  Eqn  (2-20)  and  (2-21)  into  (2-19)  to  obtain 


Atc  -  At  +  5tc  - 


where 


- -  (y  +  Vl)  +  - —  (Y2  +  v2 

1  -a2  (1  "a2  ) 


)  (2-22) 


8tc  -  (At2c  -  a2  Atlc)  /  ( 1-  a  2 ) 


Since  the  compensated  time  delay  measurement  has  been 
determined,  the  ionospheric  error  compensated  pseudo-range 
measurement  is  obtained  simply  by  multiplying  by  the  speed  of 
light,  c 


cAtc  =  cAt  +  c5tc  - 


-  {y  +  Vl)  +  - —  (y2  +  v2) 

1  -  a2  (1  -a2  ) 


-  R  +  5Rc - —  (yx  +  VX)  +  - -  (  Y2  +  V2)  (2-23) 

1  ~0T  (1  “  a2  ) 

With  this  range  measurement  model,  the  next  step  generates 
the  observation  model  for  the  range  measurement  between  the  i-th 
observer  and  the  j-th  GPS  SV. 

Let  Ps(j)  denote  the  position  vector  of  the  j-th  GPS  SV 
expressed  in  the  earth-centered  earth-fixed  (ECEF)  coordinate 
frame.  The  line-of-sight  (LOS)  vector  from  the  i-th  observer  to 
the  j-th  SV  is  then  the  vector  P^s(i,j).  This  geometry  is  shown 
for  the  two-dimensional  case  in  Figure  2.1.  From  this  figure  it 
is  seen  that 


—As ( ^ ' j )  *  Ps(3>  - 


(2-24) 


The  range  R(i,j)  is  a  nonlinear  function  of  PAs(i,j).  if  the 
computed  LOS  vector  is  defined  as  PAs,  then  the  range  is  given 
by 

R(i,j)  =  |PAs(i,j) | 

*  |PB(j)  -  PA(i) I 

=  [ ?As ( 1 ' 3 )  PAs(i,j)]1/2  (2-25) 


where  |-  |  is  the  magnitude  operator. 


FIG.  2.1.  TWO  DIMENSIONAL  NAVIGATION  MEASUREMENT  GEOMETRY 


Since  a  linearized  measurement  model  is  required  for 

A 

estimation,  R(i,j)  and  Rc(i,j)  are  expanded  in  Taylor  series 
about  the  actual  positions  to  first  order  yielding  (4:27) 


R  =  R  +  5PC  +  5  P* 

»*A 


(2-26) 


Taking  the  partial  derivatives  produces 


8R(i, j)  _  Pas ( 1 ' 3 ) 
aps(j>  R(i,j) 


—As ( 1 ' ^ ) 


3R(i, j) 

apA(j) 


R(i/ j) 


=  ~  uL(i, 


(2-27) 


(2-28) 


It  should  be  noted  that  the  first  partial  derivative  of  the 
range  is  the  unit  line  of  sight  vector.  Thus,  (2-26)  becomes 


R(i,j)  =  R(i,j)  +  UAs(i,j)  [5PS ( j )  -5PA(i)]  (2-29) 


The  linearized  measurement,  z(i,j),  can  be  formed  using 
(2-23)  and  (2-29)  as 


z(i,j)  =  Rc ( i , j )  -  R(i, j ) 


T  ca 

=  -  Has(i,3)  [5PS ( j )  -  SPA(i)]  +  SRc(i) - -  • 

1  -<*2 


C  ( i r  j )  +  v1(i,j)]  +  - -  [y2(i»J)  +  v2  ^ 1 ' ^  )  3 


1  -a2 


(2-30) 


The  errors  in  observer  position,  5PA(i),  are  the  first  three 
elements  of  the  inertial  navigation  system  error  state  model 
vectors  x2,  and  *3  as  given  by  Eqn  (2-6)  for  i=l,  2,  3 
respectively.  The  GPS  position  errors,  5Ps(j),  are  modelled  as 
unknown  biases  along  each  j-th  line  of  sight  and  are  represented 
by  the  error  state  vector  described  in  Section  2.2.2.  Thus, 
(2-30)  can  be  written 


rn  GQf 

2 ( i / j )  =  UAs(i, j)5PA(i)  +  5Rc(i) - —  [yx(i,j)  +v1(i,j)] 


1  -a 


2  L  r2 


[V?(i/j)  +  vP(i,j) ]  +  bP(j) 


(2-31) 


Note  that  the  additional  term  bE(j)  is  added  to  account  for 
unmodelled  LOS  biases. 


In  order  to  reduce  the  dimension  of  the  error  state  vector, 
the  code  loop  errors,  and  y2,  and  the  measurement  noise 
terms,  v-^  and  v2 ,  are  combined  respectively  to  obtain  two  random 
variables  Y  and  v.  By  letting 


y(i,j)  - - P  +  - 7  >2  ( i » 1 ) 

1  -  a  1  -  a2 


v(i,j)  - - ■?  v,(i,j)  +  - p  v2(i,j) 


Eqn  (2-31)  is  rewritten  as 


z(i,j)  *  UAs(i, j)5PA(i)  +  $Rc(i)  +  y(i/j)  +  v(i,j) 


b(j) 


(2 
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The  code  loop  tracking  errors  are  modelled  as  first-order 
Markov  processes  such  that  the  combined  code  loop  dynamics  are 
given  by 


1 

y  =  -  -  y  +  w 
r 


(2-33 


where  r  is  a  relatively  long  time  constant  given  by 


T  = 


(2-34 


4B 


n 


and  Bn  is  the  noise  equivalent  bandwidth  of  the  code  loop 
tracking  filter.  The  user  clock  offset  error,  5RC,  is  modelled 
as  a  constant  drift 


(2-35 


where  5fc  is  the  user  clock  frequency  error  and  5fc/fc  is  a 
random  constant  drift  bias.  Thus, 


These  two  states  model  the  clock  error.  The  equivalent  pseudo¬ 
range  bias  term,  bE,  is  modelled  as  a  random  constant  such  that 


bE  *  0  (2-37) 

With  the  navigation  measurement  model  now  complete,  the 
emitter  measurement  model  is  developed  next. 

2.4.2.  EMITTER  MEASUREMENTS 

Since  the  transmission  time  of  a  received  emitter  signal  is 
not  known  a  priori,  the  transit  time  cannot  be  used  to 
determine  the  range  to  the  emitter.  The  only  information 
available  to  each  airborne  observer  is  the  time  of  arrival  (TOA) 
of  the  emitter  signal  at  its  receiver.  For  the  case  of  the  GPS 
SVs,  the  quantities  of  interest  require  that  computations  be 
done  in  three  dimensions.  However,  for  the  case  of  the 
observers  and  the  emitter,  the  planar  projection  of  the  three 
dimensional  emitter  location  error  ellipsoid  onto  the  horizontal 
surface  is  the  quantity  of  most  significance.  This  is  because 
it  is  assumed  that  the  emitter  vertical  position  is  known  with 
sufficient  accuracy. 

Using  three  independent  measurements  of  the  emitter  signal, 
one  from  each  airborne  observer,  three  TOA's  are  obtained  with 
the  time  differences  of  arrival  ( TDOA )  defined  as: 
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tdoa12 

=  TOA  x 

-  TOA  2 

(2-38) 

tdoa23 

=  TOA  2 

-  TOA  3 

(2-39) 

The  measurement  TOA^  can  be  written  as  the  sum  of  the  actual 
TOA^  and  the  measurement  error  STOAj_ 


TOAj.  =  TOAi  +  5T0Ai 


(2-40) 


Substituting  (2-40)  into  the  form  of  (2-38) 


TDOA i  j  =  TOAi  -  TOAj  +  STOAi  -  5  TOA j  (2-41) 


Next  the  individual  error  terms  making  up  5TOA;  are  examined 
(4:31) 


5  TOA^  =  STROPi  +  SRECi  +  SCLCKi  +  vi 


(2-42] 


where 


5TROP^  =  time  delay  due  to  errors  in  tropospheric  model 


fiREC.; 


5 CL CKi 


time  delay  due  to  uncalibrated  receiver  delay 
observer  clock  bias  error 
measurement  corruption  noise 


Next,  Equation  (2-42)  is  substituted  into  (2-41)  yielding 


TDOAij 


(2-43 


=  TOA^  -  TOAj  +  SCLCKj^j  +  STROPj^  -  5  TROP j 
+  5 REC^  -  SRECj  +  vj_j 


where 


and 


5  CLCK^  j  =  SCLCK  j_  -  SCLCKj 


v 


ij 


=  Vj 


V 


j 


The  time  difference  of  arrival  in  (2-43)  is  scaled  by  the  speed 
of  light  to  yield  the  range  difference 


AR(i,j)  -AR(i,j)  +  8Rc(i)  -  5Rc(j)  +S RT(i)  -6RT(j) 

+  SRR(i)  -  5RR(j)  +  v(i,j)  (2-44 


where 


AR(i,j) 

= 

cfTOAi  -  TOAj) 

AR(i, j) 

= 

cfTOAi  -  TOAj) 

5Rc(i) 

= 

c  8 CLCK^ 

5  RT ( i ) 

= 

C  5  TROPj^ 

5RR(i) 

= 

c  §REC^ 

v(i, j) 

c(wi  -  Vj) 

It  should  be  noted  that  each  is  modelled  as  a  white  Gaussian 
noise  of  mean  zero  and,  further,  that  the  V; 's  are  independent 


of  each  other  such  that 


Vi  =  N  [ 0 , R] 


v(i,j)  =  N  [0,2R] 


(2-45) 


where  R  is  the  variance  of  each  Vi  (4:32) 


Note  that  the  clock  is  assumed  to  be  the  same  for  the  entire 
emitter  location  problem.  Therefore,  the  same  clock  error 
states  are  shared  for  the  navigation  as  well  as  the  emitter 
solution.  The  error  terms  of  5RC,  5Rt,  and  5rr  are  related  to 
the  error  state  vector  x(t)  by  noting  that  5RC  is  represented  by 
cx5  of  Section  2.2.2,  5Rt  is  represented  by  cxi;L  and  5rr  is 
represented  by  cx10  with  x10  and  xn  described  in  Section  2.3. 


As  with  the  GPS  measurement  development  in  Section  2.4.1, 
the  next  step  develops  the  observation  model  from  the  range 
measurement  model. 


The  error  state  observation  is  well  modelled  by 


where 


z(i,  j)  =  AR(i,  j)  -  AR(i,  j) 


AR(i,j)  =  computed  or  estimated  measurement 


(2-46) 


V-  .'.V- .V ■  .V 


Let  PE  denote  the  position  vector  of  the  emitter  expressed  in 

a  ,  A 

the  ECEF  coordinate  frame.  Next,  let  PA(i)  and  represent 

the  position  vectors  of  the  i-th  and  j-th  observer  aircraft 
respectively,  again  expressed  in  the  ECEF  frame.  The  LOS 
vectors  from  the  emitter  to  the  i-th  and  j-th  observer  are  then 

A  A 

the  vectors  PAE(i)  and  — AE ( i )  /  respectively.  This  measurement 
geometry  is  described  in  two  dimensions  in  Fig.  2.2.  It  is 
noted  from  the  figure  that 


PaE^)  =  2e  -  {2-Al 

A 

So  the  range  difference  AR(i,j)  is  a  nonlinear  function  of 

A  A 

PAEU)  and  PAE(j).  Using  these  computed  lines  of  sight  vectors 
allows  the  range  difference  to  be  expressed  as 


AR(i,j) 


IPae^)!  '  ISae<3)  I 

A  A  A 

|PE  -  PAU)  I  -  |PE  - 


PA(i)l 


(2-48 


A  A 

Note  that  PA(i)  and  PA(j)  are  those  estimated  quantities  which 
have  been  used  in  Equation  (2-25)  for  deriving  the  navigation 
measurement  model.  Again,  a  linearized  measurement  model  is 

A 

required  so  AR(i,j)  is  expanded  in  Taylor  series  to  first  order 


about  the  actual  observer  and  emitter  positions  yielding 


A  3ART(i,j)  5AR1  (i,  j ) 

AR(i,j)  =  AR(i,  j )  +  -  5PE  +  -  5 P A ( i ) 

9PE  SPA(i) 


3ART(i,  j) 

+  -  SpA(j)  (2-49) 

3PA(j) 


Evaluating  each  of  the  required  partial  derivatives  using  (2-48) 
yields 


3ART(i,j)  3 

-  -  -  [  IEaeCMI  ~  I  — AE C ^ )  I  ]  (2-50) 

9PE  spE 


Therefore,  the  partial  derivatives  of  (2-49)  are  written  as 

3ART(i,  j)  PAE(i)  —  AE ( j  ) 

dPE  Pj 

-  UAE^)  -  —  AE  ( 3  )  (2“ 

Again,  note  that  the  derivatives  of  the  magnitudes  of  the  LOS 
vectors  conveniently  become  the  unit  line  of  sight  vectors;  i 
this  case  from  the  i-th  and  j-th  observers  to  the  emitter. 

The  third  term  on  the  right  hand  side  of  (2-49)  is  evaluated 

3AR(i,j)  _  PAE(i) 

SPA(i)  Rj. 
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Finally,  the  last  term  in  (2-49)  is 


3AR(i,j)  _  PAE(j) 

SPA(j)  Ri 

=  UATE(j)  (2-53) 

Substituting  (2-44)  and  (2-49),  with  the  partial  derivatives 
evaluated,  into  (2-46)  produces  the  required  error  measurement 
equation 


z(i,j)  -  5Rc(i)  -  SRc(j)  +  5RT(i)  -  5 RT ( j ) 

+  8 Rr ( i )  -  5RR(j)  +  v(i,j) 

-  [UAE(i)  -  UAE(j)]  $PE  +  UAE(i)5PA(i) 

-  UAE(j)5PA(j)  (2-54) 

This  completes  the  model  development  for  system  error 
dynamics,  including  the  INS,  GPS  and  emitter  location  system 
errors,  as  well  as  the  measurement  models  for  navigation  and 
emitter  measurements.  In  the  next  chapter  the  estimator,  cost 
function,  and  cost  gradient  are  developed. 
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III.  ESTIMATOR  AND  COST  FUNCTION  DEVELOPMENT 


3.1  INTRODUCTION 


The  models  developed  in  Chapter  2  are  the  cornerstones  upon 
which  an  estimator  can  be  built.  The  purpose  of  this  estimator 
is  to  provide  a  time  history  of  the  statistical  properties  of 
the  error  state  vector  x(t)  as  it  is  propagated  forward  in  time 
from  some  initial  condition  x(tQ)  =  xQ.  For  this  study,  a 
discrete-time,  time-varying  estimator  is  used.  Because  the 
measurement  model  is  nonlinear,  an  extended  Kalman  filter  was 
selected  and  is  described  in  Section  3.2  of  this  chapter. 


In  order  to  develop  a  cost  function  for  an  emitter  locating 
system,  it  is  necessary  to  examine  the  figures  of  merit 
available  for  evaluating  the  performance  of  such  a  system.  The 
most  obvious,  and  most  crucial,  performance  criterion  is  the 
accuracy  with  which  the  system  can  estimate  the  actual  emitter 
location.  Although  the  accuracy  of  this  estimate  can  be 
expressed  in  any  of  several  ways,  the  mean  squared  miss  distance 
serves  as  the  basis  for  cost  function  development  in  this  study. 
The  development  of  the  cost  function  and  some  of  its  limitations 
are  described  in  Section  3.3.  Using  the  derived  cost  function, 
the  cost  gradient  is  developed  in  Section  3.4.  Finally,  the 
minimum  cost  search  technique  used  in  this  study  is  presented  in 
Section  3.5. 


ESTIMATOR  DEVELOPMENT 


In  order  to  use  the  Kalman  filter  as  the  estimator  for  this 
study,  certain  assumptions  need  to  be  made.  These  assumptions 
are  that  (1)  the  state  dynamics  model  and  the  measurement  model 
are  linear,  (2)  system  dynamic  driving  noise  and  measurement 
corruption  noise  processes  are  well  described  by  Gaussian 
processes,  and  (3)  all  system  noise  processes  are  well  modelled 
as  "white".  With  these  assumptions,  the  extended  Kalman  filter 
is  developed  (6) . 


First,  the  nonlinear  error  state  dynamics  model  given  by 
(2-2)  is  expanded  in  Taylor  series  about  the  current  best 
estimate  of  the  error  state  vector  and  truncated  to  first  order 
terms 


A  df (X, t) 

f(X,t)  =  f[x(t),t)  +  - 


(3-D 


d X  I  X  =  x(t) 


where  8x  represents  the  error  state  vector.  Defining 


df(x,t)  | 

F(t)  =  -  | 

dx  |  x  =  x  (t ) 


(3-2) 


produces  the  linearized  dynamics  model 


x  =  F  (t )  x  +  G  w 


(3-3) 


Note  that  for  the  first  time  interval,  the  nominal  point  for 

.  A  A 

expansion  is  x(tQ)  =  xQ.  Also  note  that  the  model  is 
relinearized  about  the  new  estimate  each  time  it  is  computed. 
This  allows  the  nominal  trajectory  to  be  updated  on  a  continuing 
basis  to  ensure  that  deviations  from  the  nominal  remain  small. 

Using  the  same  approach  for  h(x,t),  the  nonlinear  measurement 
vector  given  by  (2-4)  is  linearized.  Define 


ah(x,t)  I 

H  ( t)  =  -  |  A  (3-4) 

ax  |x  =  x(t) 

resulting  in  the  linearized  measurement  equation 


zfti)  =  H(ti)x  +  v  (3-5) 

Forming  the  equivalent  discrete-time  system  model,  Eqn  (3-3) 
becomes  the  stochastic  difference  equation 


x^i+l)  “  4>(ti+1,ti)  xfti)  +  w^t^ 


(3-6) 


where 


E  [wd(ti) ]  =  0 

E  [w(j(ti)  wd(ti)]  -  Qd^ti+l^i) 


(3-7) 


and 


X 


Qd(t,ti)  =  <t(  t,r)  G(r)  Q(r)  GT  (  r  )  <$T  (t ,  r  )  dr 


(3-8) 


Note  that  <t>(t,t±)  is  the  state  transition  matrix  such  that 


$(t0,t0)  =  I  (3-9) 

$>(t2,t0)  =  c^(t2 , tx)  <Mt1,t0)  (3-10 

The  measurement  equation  (3-5)  already  models  discrete-time 
measurements,  so  it  remains  unchanged. 

Having  satisfied  the  three  initial  assumptions  in  this 
section,  the  Kalman  filter  equations  are  stated. 

The  estimate  is  propagated  forward  in  time  using  (6) : 

A  -  A  +  a 

x(ti+1)  «  X(tA)  +  /  f[x(t/ti),t]  dt 

1 

P(t-+1)  =  <D(ti+1,ti)  P(tJ)  <t>T(ti+1,ti)  +  Qd(ti)  (3-11 

where 

P(ti  =  0)  =  PQ  (3-12 

and  the  measurement  update  equations  are  given  by  (6) : 

K(tA)  =  P(tJ)  H1^)  [H(ti)  P(t-)  HT(ti)  +  R(ti)]"1  (3-13 

X(t-)  =  X(tJ)  +  K(ti)  [Zi  -  H(tA)  X(t-)]  (3-14 

P(ti)  =  P(t l)  -  K(ti)  H(ti)  P(t-) 


( 


Using  the  assumptions  presented  at  the  beginning  of  this 
section  allows  the  conditional  probability  density  function, 


fx/z(<£,Z)  ,  to  be  completely  described  by  only  the  first  two 
moments.  This  provides  mathematical  tractability  and  allows  a 
time  history  of  x(t)  to  be  generated,  conditioned  on  all 
measurements  up  through  time  t,  by  propagating  only  the 
conditional  mean  and  the  covariance.  The  second  moment,  the 
error  covariance  matrix  P(t) ,  is  used  for  the  error  covariance 
analysis  in  this  study. 

3.3  COST  FUNCTION  FORMULATION 


In  general,  the  figure  of  merit  used  to  evaluate  the 
performance  of  an  emitter  locating  system  is  the  accuracy  with 
which  the  system  can  estimate  the  actual  emitter  location  in 
three  dimensions.  For  this  study,  the  problem  reduces  to  an 
essentially  two-dimensional,  or  planar  problem,  due  to  the 
assumption  that  the  vertical  emitter  position  is  known 
precisely.  Therefore,  the  quantity  of  real  interest  is  the 
emitter  location  error  ellipse  formed  by  the  projection  of  the 
error  ellipsoid  onto  the  x-y  plane.  The  error  ellipse  used  in 
this  study  is  the  circular  error  probable  (CEP)  contour,  which 
is  a  contour  of  equal  joint  probability  density  that  yields  a 
probability  of  0.5  when  integration  of  the  probability  density 
function  is  carried  out  over  the  area  enclosed  by  the  contour. 


The  expression  used  for  the  CEP  in  this  study  is  (4:40) 


CEP  =  0.5887  (£7X  +  oy)  (3-16) 

Note  that  this  expression  represents  the  CEP  for  an  actual 
circle  where  ax  and  ay,  the  semi-major  and  semi-minor  axes,  are 
equal.  For  the  more  realistic  case  where  the  equal  probability 
density  contour  is  an  ellipse  with  ax  ?  ay,  it  is  found  that 
(3-16)  approximates  the  true  CEP  to  within  three  percent  as  long 
as 


0.15  <  —  <  1.0  (3-17) 

a 

y 

Because  of  the  validity  of  the  approximation  given  by  Eqn 
(3-16) ,  the  cost  function  may  be  defined  in  terms  of  CEP.  Since 
the  CEP  is  a  function  of  ox  and  ay,  the  error  covariance  matrix, 
P(t)  ,  is  examined  to  see  how  these  quantities  can  be  obtained. 
First,  the  two  states  of  interest  are  the  emitter  location 
errors  in  the  horizontal  plane,  x(83)  and  x(84).  The  four 
elements  of  P(t)  needed  to  obtain  ctx  and  cry  are  given  by 

I  Pt (83,83)  Pt (83,84)  | 

pemit  =  I  I  ( 3  —  1 8 ) 

|  Pt(84,83)  Pt(84,84)  | 


or  in  simplified  notation 


pemit 


pll 

P12 

P2 1 

P22 

(3-19) 


Note  that  the  eigenvalues  of  Pem^t  are  the  magnitudes  of  the 
principal  axes  of  the  error  ellipse  and  the  eigenvectors 
represent  the  principal  axes  directions.  If  it  is  assumed  that 
the  principal  axes  are  rotated  only  a  very  small  amount  from  the 
coordinate  axes,  then  the  cross  covariance  terms  P12  and  P21  are 
approximately  zero  and  the  diagonal  terms  are  given  as 


^11  —  *7^  and  ?22  —  °y 


(3-20) 


Now  it  is  possible  to  define  a  mean  squared  miss  distance  cost 
function,  J,  as 


J(PS)  =  ax  +  °y 


~  P11  +  P2  2 


(3-21) 


which  is  related  to  the  root  mean  square  (RMS)  miss  distance  by 

RMS  =  J1/2  (3-22) 

An  alternate  method  is  to  define  the  cost  function  in  terms  of 
the  CEP  such  that  JCEp  is  given  by 


'<<:  >j|  ■j.'W  n  >>>  9  %  r«l  rv* 


a 


JCEP^— S)  =  ax  +  ay  (3-23) 

which  is  directly  related  to  CEP  by 

CEP  =  0.5887  JCEp(Ps)  (3-24) 

As  pointed  out  by  Lewantowicz  (4),  the  cost  gradient  for  mean 
squared  miss  distance  computation  g(Ps),  is  easier  than  the 
computation  for  the  CEP  cost  gradient,  ScEP^-s^*  For  th^s 
reason,  J(PS)  given  by  (3-21),  the  mean  squared  miss  distance 
cost,  is  the  cost  function  which  is  actually  minimized  in  this 
study.  The  results  are  presented  in  terms  of  both  CEP  and  RMS 
miss  distance. 


Note  that  the  error  covariance  matrix,  P(t^) ,  is  iterated 
from  the  initial  conditions  to  steady-state  values  at  each  step 
of  the  cost  search.  This  prevents  the  accumulation  of  past  SV 
position  related  information  in  the  cross-covariance  terms  of 
the  P  matrix  as  the  satellite  geometry  changes. 


Using  the  cost  function  described  by  (3-21),  the  cost 
gradient  is  developed. 


(3-25) 


J(PS)  -  tr  [  E  Pt  ] 

=  Pt (83,83)  +  Pfc (84,84) 


where 

E  is  a  square  91  x  91  matrix  of  all  zeros,  except  for  l's 
in  the  83rd  and  84th  diagonal  positions. 

and 

Pt  is  the  error  covariance  matrix  at  time  t 

From  (3-25)  and  (3-11),  it  can  be  seen  that  the  cost  J  is  a 
function  of  the  vector  Ps  described  in  Section  2.4.1.  This  Pg 
starts  out  as  a  12  dimensional  vector  describing  the  three- 
dimensional  position  of  each  of  4  satellites  in  the  ECEF  frame. 
Note  that  the  vector  Ps  is  constrained  by  the  fixed  orbital 
radius  of  the  GPS  satellites.  For  the  j-th  satellite,  the 
vector  is  constrained  such  that  (4:68) 

Xj  =  (Rg  -  Vj  -  Zj  )1/2  (3-26) 

The  geometry  chosen  is  such  that  the  x-axis  of  the  ECEF  frame 
passes  through  the  "center"  of  the  observer  geometry,  where  the 
center  is  defined  as  a  point  equidistant  from  the  endpoints  of  a 
line  segment  connecting  observer  1  and  observer  3.  The  y-z 
plane  forms  the  two-dimensional  space  onto  which  the  SV  and 
observer  positions  are  projected.  Note  that  for  projection 
purposes,  the  center  of  the  observer  geometry  is  located  along 
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the  x-axis  at  an  altitude  h  above  the  surface  of  the  earth. 
Incorporating  the  effect  of  this  constraint  yields  an  eight 
dimensional  parameter  vector  Ps. 

— s  =  [  Yl'  zl>  Y2>  z2 '  Y2>  z3 /  y4'  2 4  ]  (3-27) 

For  each  SV,  this  two-dimensional  space  is  further 
constrained  by  the  initial  assumption  that  satellites  used  must 
lie  at  or  above  five  degrees  elevation  from  the  local  horizon 
measured  from  the  intersection  of  the  y-z  axes.  The  result  of 
this  constraint  maps  all  allowable  satellite  positions  into  a 
circle  of  radius  Rc  centered  at  the  origin  of  the  y-z  coordinate 
system. 


Fig.  3.1  depicts  a  triangle  with  angles  A,  B,  and  C  with  sides 
opposite  these  angles  given  by  a,  b,  and  c  respectively.  From 
this  figure  it  can  be  seen  that 


R_  |  =  R_  cos  /3 


(3-28; 


Referring  to  Fig.  3.1,  find  the  angle  C  using  the  law  of  sines 


Rg  Rg+h 


sin  B  sin  C 


(3-29) 


where  B  =  95°. 


solving  for  the  angle  C  yields: 


1  «*t  t  l  I 


I»>  .  .»>  »u- A  t«fc 


c 


(3-30) 


=  arcsin 


(R.+h)  sin  95c 


Now  find  /3 


(3  =  90°  -  A 


(3-31) 


where 


A  =  180°  -  (B  +  C) 

Substituting  Rs  =  8.4116  x  107  into  Eqn  (3-28)  yields 


(3-32) 


|RC!  =  7.934  X  107  feet 


(3-33) 


This  figure  is  now  used  to  determine  when  projected  satellite 
positions  violate  the  LOS  elevation  constraint.  If  the 
satellite  position  projection  violates  this  constraint,  then  it 
is  reset  onto  the  constraint  boundary. 


The  actual  cost  gradient,  g(Ps),  is  developed  in  Appendix 
A.  This  cost  gradient  is  now  incorporated  into  a  minimum  cost 
search  as  described  in  the  next  section. 


3.5  MINIMUM  COST  SEARCH  TECHNIQUE 

The  search  technique  selected  for  minimizing  the  cost 
function  in  this  study  is  a  steepest  descent  algorithm.  The 
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search  scheme  is  based  on  the  premise  that  the  computed  gradient 
g(Ps)  indicates  the  direction  of  maximum  rate  of  increase  in  the 
cost  for  a  given  Ps.  Therefore,  the  negative  gradient  direction 
must  represent  the  direction  of  maximum  rate  of  decrease  in  the 
cost.  It  is  possible  to  proceed  along  the  negative  gradient 
direction  until  a  minimum  cost  along  this  gradient  is 
determined.  At  this  point  a  new  gradient  direction  is  computed 
and  search  proceeds  along  the  new  negative  gradient  direction. 
The  search  continues  until  the  gradient  vector  computed  has  a 
magnitude  smaller  than  a  preset  threshold.  The  initial  step 
size  for  the  search  is  chosen  by  trial  and  error.  It  is 
important  to  consider  the  trade-offs  between  the  number  of 
intermediate  points  reguired  to  reach  the  minimum  along  a 
gradient  direction  and  the  possibility  of  grossly  overshooting 
that  minimum  when  determining  the  initial  step  size. 

The  minimum  cost  function  point  along  a  given  gradient 
search  direction  gQ  is  found  by  monitoring  the  angle  condition 
between  gQ  and  g^  using  the  cosine  of  the  angle  between  the  two 
vectors.  This  cosine  is  computed  by  using  the  definition  of  the 
inner  product  (or  vector  projection) 

<gQ  /  ?i>  =  I  g0  l  l^i  I  cos  e  (3-34) 

where  6  is  the  angle  between  the  two  vectors. 
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Solving  for  cos  6 


cos  0  = 


<So  '  Ii> 


(3-35) 


The  minimum  cost  point  along  gQ  occurs  where  g^  is  orthogonal  to 
gQ,  i.e.  where  cos  6  =  0.  When  this  orthogonality  condition  is 
met,  the  g^  becomes  the  new  gQ  and  Ps^  becomes  PSo  as  shown  in 
Fig.  3.2. 

The  vector  Ps^  is  computed  as 


-Si  =  £so  +  si  1c 


where 


Sj_  =  (  1  +  0.9  cos 


(3-36) 


(3-37) 


Note  that  SQ  is  the  initial  step  size  selected  and  that  the 
algorithm  automatically  reduces  the  step  size  as  the  gQ,  g^ 
orthogonality  condition  is  approached. 

This  search  algorithm  is  modified  due  to  the  fact  that  the 
satellites  may  move  in  either  an  unconstrained  or  constrained 
manner.  Constrained  movement  is  required  when  the  gradient 
search  direction  indicates  the  satellite  should  be  moved  onto  or 
beyond  the  Rc  boundary.  In  this  study,  the  constrained  and 
unconstrained  satellite  movements  are  computed  independently 
using  separate  initial  step  size  parameters.  Now,  gQ  is  given 
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(3-38) 


<?Q  =  t^CON  +  [9Q] 


oJ UNCON 


Whenever  a  satellite  reaches  the  Rc  boundary,  gQ  is  redefined 
and  the  search  begins  along  the  new  gradient  search  direction. 
The  [gD] uncon  redefined  using  g^  for  the  satellites 

that  have  not  reached  the  boundary  and,  similarly,  Ps^  becomes 
PSo.  The  [go]C0N  is  actually  the  projection  of  g^,  for  the 
satellites  on  the  boundary,  onto  a  vector  tangent  to  the 
boundary  circle  at  Ps^.  This  projection  is  done  so  that,  in 
approximation,  gradient  search  is  done  along  the  circular 
boundary  as  indicated  in  Fig.  3.3.  If,  during  gradient  search 
along  the  boundary,  g^  indicates  satellite  movement  should  be 
off  the  boundary,  then  gQ  is  again  reset  such  that  the 
appropriate  satellite  may  move  off  the  boundary  constraint. 


To  find  [go]C0N  first  define  the  normal  vector  at  Ps^  for 
the  j-th  satellite 


n  -  [Yj  »  Zj] 


(3-39) 


The  tangent  vector  at  Ps^  is  then 


t  =  [Zj  ,  -yj] 


which  is  normalized  to  form  the  unit  tangent  vector 


(3-40) 
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The  projection  angle  0  is  given  by 


cos  0 


<2ij 


ut; 


lij 


(3-42) 


The  projection  magnitude  is  computed  as 

I  t^olcON  I  “  I  Sfij  I  cos  0  (3-43) 

Finally,  gG  for  the  constrained  satellite  is  given  as 

[Sfo^CON  =  *  ^ol  CON  I  -t  (3-44) 

The  termination  condition  for  search  in  the  case  of  constrained 
satellites  can  obviously  no  longer  be  based  on  the  magnitude  of 
the  computed  gradient  vector  being  small.  The  appropriate 
termination  condition  in  this  case  is  when  the  computed  gradient 
vector  is  orthogonal  to  the  tangent  vector  or,  in  other  words, 
when  the  projection  of  the  gradient  vector  oijto  the  tangent  is 
very  nearly  zero. 

Having  developed  the  basic  mechanisms  for  estimation,  cost 
function  computation,  and  minimum  cost  search,  the  analysis  is 
accomplished.  The  next  chapter  reports  the  results. 


IV.  RESULTS 
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4.1  INTRODUCTION 


In  this  study,  the  work  accomplished  by  Lewantowicz  in  his 
December  1985  Master  of  Science  thesis  (4)  is  expanded  and  used 
as  a  starting  point  for  further  research.  In  particular,  the 
effect  of  iterating  the  error  covariance  matrix  to  steady  state 
after  each  satellite  repositioning  movement  is  examined.  This 
chapter  examines  the  optimum  positions  of  GPS  satellits  for  the 
four,  three,  and  two  satellite  cases  such  that  minimum  errors  in 
the  estimate  of  emitter  location  are  made.  The  results  of  the 
four  and  three  satellite  cases  are  compared  directly  with  the 
results  obtained  by  Lewantowicz.  In  addition,  the  validity  of 
using  a  fourth  satellite  directly  overhead  to  simulate  a  three 
satellite  case  where  each  observer  has  a  precise  clock  is 
examined.  The  performance  of  the  emitter  location  system  is 
evaluated  for  each  optimum  satellite  configuration  and  the 
results  are  compared. 

Using  the  satellite  geometry  relations  to  the  cost  data 
obtained  in  all  the  optimization  runs,  the  basic  guidelines  are 
developed  for  a  satellite  selection  algorithm  formulation.  The 
factors  considered  in  forming  the  algorithm,  as  well  as 
justifying  data,  are  presented  in  section  4.5. 
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FOUR  SATELLITE  OPTIMIZATION 


The  four  satellite  position  optimization  problem  begins 
with  the  SVs  and  observers  in  the  initial  configuration  shown  in 
Fig.  4.1.  These  positions  were  selected  to  correspond  with 
those  used  by  Lewantowicz  so  that  a  direct  comparison  can  be 
made.  Note  not  only  the  initial  geometry  of  the  problem,  but 
also  the  characteristics  of  each  of  the  error  ellipses.  Each 
observer  error  ellipse  has  a  particular  size  and  orientation, 
and  the  emitter  error  ellipse  has  important  characteristics  of 
size  and  elongation,  which  are  ultimately  reflected  in  the  CEP. 
It  should  be  pointed  out  that  the  elongation  of  the  emitter 
error  ellipse,  in  general,  is  not  taken  into  account  if  the 
chosen  satellite  geometry  only  minimizes  the  observer  navigation 
problem.  For  many  practical  applications,  the  orientation  and 
elongation  of  the  emitter  error  ellipse  is  of  critical 
importance,  especially  in  the  case  of  highly  elongated  ellipses. 
The  initial  satellite  position  data  is  given  in  Table  4.1  and 
initial  error  ellipse  parameters  are  presented  in  Table  4.2. 

From  this  starting  point,  the  goal  is  to  accomplish  a 
gradient  search  which  will  move  the  satellites  to  an  orientation 
yielding  minimum  emitter  MSMD  (and  CEP) .  However,  before  each 
of  the  gradient  computations  is  performed,  it  is  desirable  to 
iterate  the  error  covariance  matrix  to  steady-state  at  a 
constant  SV  position  to  examine  the  Kalman  filter  convergence 
time  history. 


Table  4.1.  Four  Satellite  Initial  Positions 


1 

I 

Azimuth 

Elevation  (— ) 

LOS  Bias 

SV-1 

1 

1 

1 

35.0 

69.0 

9.8 

SV-2 

1 

1 

228.0 

31.9 

9.6 

SV-3 

1 

i 

1 

86.0 

18.6 

10.7 

SV-4 

1 

1 

24.6 

12.2 

9.7 

Table  4.2.  Four  Satellite  Initial  Error  Ellipse  Parameters 

1  Axis  1  ( ft)  Axis  2  Angle  from  Y  X~1  CEP  ( f 


Although  many  methods  exist  to  determine  the  number  of 


iterations  required  to  reach  steady-state  operation,  this  study 
examined  the  relationship  between  the  number  of  iterations  of 
the  error  covariance  matrix  and  the  parameters  describing  the 
emitter  error  ellipse.  In  particular,  the  change  in  magnitude 
of  CEP,  ACEP,  from  iteration  i  to  i+1  is  computed  and  the 
minimum  value  is  sought  as  a  function  of  i.  In  addition,  the 
convergence  angle,  which  is  the  angle  between  the  gradient 
vector  computed  at  iteration  i  and  the  gradient  vector  at 
iteration  i+1,  is  monitored  to  determine  when  the  computed 
gradient  vectors  become  coaligned.  The  ACEP  for  200  Kalman 
filter  iterations  of  the  error  covariance  matrix  propagation  and 
update  is  presented  in  Fig.  4.2.  Note  the  rapid  drop  in  CEP 
improvement  as  the  iterations  increase,  with  a  local  minimum 
occuring  at  eight  iterations.  This  minimum  is  clearly  shown  on 
the  expanded  iteration  scale  in  Fig.  4.3.  The  behavior  of  ACEP 
after  iteration  eight  is  presented  on  an  expanded  ACEP  scale  in 
Fig.  4.4.  Note  the  transient  dynamic  behavior  of  the  matrix 
Ricatti  equation  solution  as  reflected  by  this  CEP  measure  in 
Fig.  4.4,  and  the  second  local  minimum  at  approximately  35 
iterations.  From  the  data  presented,  the  obvious  choice  for  the 
number  of  iterations  of  the  error  covariance  matrix  to  use  for 
steady  state  would  be  one  of  the  two  local  minima.  Due  to  the 
computational  loading  required  to  propagate  and  update  the  error 
covariance  matrix  at  each  cost  function  search  point,  and  the 
magnitude  of  improvement  expected  with  increased  iterations,  the 
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Fig.  4.4.  CEP  Transient  Behavior 
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local  minimum  of  eight  iterations  was  selected  to  approximate 
steady-state  operation  for  this  study.  Gradient  vector 
convergence  at  8  iterations  is  good,  with  a  convergence  angle  of 
less  than  0.6  degrees.  With  eight  iterations  established  as  the 
number  of  iterations  prior  to  each  gradient  vector  computation, 
the  minimum  cost  search  for  the  four  satellite  case  is 
conducted.  Note  that  the  eight  iteration  figure  was  obtained 
with  the  satellites  and  observers  only  at  their  initial 
positions.  Therefore,  this  figure  may  or  may  not  be  appropriate 
for  all  possible  scenarios. 

The  satellite  positions  determined  by  the  minimum  cost 
search  are  presented  graphically  in  Fig.  4.5,  with  position  data 
and  error  ellipse  parameters  given  in  Tables  4.3  and  4.4, 
respectively.  The  paths  the  satellites  "traveled"  during  the 
gradient  search  in  reaching  their  final  positions  are  shown  in 
Fig.  4.6.  The  minimum  emitter  CEP  found  is  41.22  feet,  down 
from  the  initial  CEP  of  58.0  feet.  Comparing  Tables  4.2  and 
4.4,  it  can  be  seen  that  the  observer  horizontal  navigation 
errors  have  been  significantly  reduced.  This  reduction  results 
from  a  marked  decrease  in  LOS  bias  errors  as  indicated  in  Table 
4.3.  The  reduction  in  these  east-west  observer  navigation 
errors  directly  affect  the  shape  and  orientation  of  the  emitter 
ellipse.  In  fact,  these  reductions  are  the  driving  force  behind 
the  reduction  in  the  semi-major  axis  of  the  emitter  ellipse  and, 
ultimately,  the  emitter  CEP.  Note  that  three  of  the  four 


Table  4.3 


Four  Satellite  Final  Positions 


1 

1 

Azimuth  (-) 

Elevation  (-) 

LOS  Bias 

SV-1 

1 

1 

| 

275.3 

25.0 

2.8 

SV-2 

1 

1 

1 

277.4 

24.3 

2.9 

SV-3 

1 

1 

1 

101.5 

5.0 

2 . 6 

SV-4 

1 

1 

13.3 

5.0 

1.8 

Table  4.4.  Four  Satellite  Final  Error  Ellipse  Parameters 


1 

1 

Axis  1 

(ft)  Axis  2 

Ancile  from  Y  (2) 

CEP  (ft) 

Obs-1 

1 

1 

1 

3.15 

2.02 

-41.5 

3.04 

Obs-2 

1 

1 

1 

3.16 

2 . 02 

-42.0 

3.05 

Obs-3 

1 

1 

1 

3.00 

1.77 

-35.5 

2 .81 

Emit 

1 

1 

62.98 

7.05 

-  0.07 

41.22 

♦TMtW 


(  FT 


satellites  tended  to  align  themselves  along  the  east-west  axis 
in  opposition  to  each  other.  This  alignment  produced  the  noted 
reductions  in  the  east-west  emitter  location  error,  and  reflects 
the  type  of  optimum  geometry  expected  at  the  outset  of  this 
study.  The  fact  that  the  semi-major  axis  of  the  emitter  error 
ellipse  lies  along  the  east-west  axis  is  a  function  of  the 
observer  geometry. 

Comparing  the  results  of  the  four  satellite  optimization 
with  those  obtained  by  Lewantowicz  (4)  indicates  that  the 
satellite  positions  are  within  10°  of  each  other  in  azimuth  and 
elevation  and  the  minimum  emitter  CEPs  are  within  two  feet  of 
one  another.  This  supports  the  Lewantowicz  hypothesis  described 
in  Section  1.5.  Next,  the  three  satellite  optimization  is 
performed. 

4.3  THREE  SATELLITE  OPTIMIZATION 

As  pointed  out  in  Section  3.3,  the  cost  function  used  in 
this  study  was  formulated  for  an  essentially  two-dimensional 
problem.  This  means  that  the  problem  is  uniquely  solved  with 
only  three  satellites  and  the  fourth  satellite  only  provides 
additional  horizontal  emitter  position  information.  The 
performance  of  the  emitter  locating  system  is  evaluated  for  two 
separate  three  satellite  cases;  first,  one  satellite  is  fixed 
overhead  and  three  satellites  are  allowed  to  move  (see  Section 


1.5)  and,  second,  the  measurements  from  only  three  satellites 
are  incorporated  into  the  error  covariance  computation. 

The  "pseudo"  three  satellite  case,  as  proposed  by 
Lewantowicz  (4)  uses  four  satellites  with  one  fixed  directly 
overhead  and  is  intended  to  simulate  the  three  satellite  case 
when  each  observer  has  a  precise  clock.  The  results  of  the 
pseudo-three  satellite  optimization  are  presented  in  Fig.  4.7 
and  Tables  4.5  and  4.6.  The  paths  the  satellites  "traveled"  in 
reaching  their  final  positions  are  shown  in  Fig.  4.8.  Note  the 
final  CEP  achieved  in  the  pseudo-three  satellite  case  is  40.71 
feet  compared  to  41.22  feet  obtained  in  the  four  satellite  case. 
These  results  seem  unreasonable  since  they  indicate  that  the  CEP 
improved  when  less  position  information  was  available.  What  has 
happened  is  that  a  three  satellite  optimization  has  not  been 
performed,  but,  rather,  a  constrained  four  satellite  case  has 
been  optimized.  The  gradient  search  performed  for  the 
constrained  satellite  case  found  a  different  local  minimum  than 
that  found  in  the  original  four  satellite  case.  The  existence 
of  this  second  local  minimum  indicates  that  finding  a  global 
minimum  using  the  gradient  search  algorithm  of  this  study  is 
highly  unlikely.  However,  the  magnitude  of  the  second  local 
minimum  varies  only  slightly  from  the  first  minimum,  indicating 
that  it  is  reasonable  to  expect  other  local  minima  to  result  in 
nearly  the  same  costs.  Note  that  satellites  1  and  2  in  the 
pseudo-three  satellite  case  have  assumed  approximately  the  same 


Table  4.5.  Pseudo  Three  Satellite  Final  Positions 


1 

1 

Azimuth  (—) 

Elevation  (-) 

LOS  Bias 

SV-1 

1 

1 

1 

277.4 

16.7 

1.6 

SV-2 

1 

1 

1 

272.3 

17.1 

1.7 

SV-3 

1 

1 

47.7 

5 . 0 

1.6 

Table  4.6 


Pseudo  Three  Satellite  Error  Ellipse  Parameters 


1 

1 

Axis  1 

(ft)  Axis  2 

Anale  from  Y  (-) 

Obs-1 

1 

1 

1 

2.46 

1.55 

-88. 0 

Obs-2 

1 

1 

1 

2.50 

1.54 

-87 . 7 

Obs-3 

1 

1 

1 

1.74 

1.35 

65.1 

Emit 

1 

1 

62.33 

6.82 

-  0.01 

CEP  (ft 
2.34 
2 . 38 
1 . 82 
40.71 


positions  that  they  did  in  the  four  satellite  case,  adding 
credence  to  the  idea  of  a  constrained  four  satellite  case. 

The  fact  that  the  pseudo-three  satellite  optimization 
yields  better  results  than  the  four  satellite  case  indicates 
that  the  simulation  of  a  three  satellite  case  using  four 
satellites,  with  one  satellite  fixed  overhead,  is  not 
appropriate.  This  is  further  verified  by  the  following  true 
three  satellite  optimization. 

For  the  true  three  satellite  case,  the  measurements  from 
only  three  satellites  are  incorporated  in  the  error  covariance 
matrix.  The  three  satellite  optimization  results  are  given  in 
Fig.  4.9  and  Tables  4.7  and  4.8.  Again,  satellite  tracks  are 
indicated  in  Fig.  4.10. 

The  results  in  Table  4.8  are  reasonable  for  the  three 
satellite  case,  with  the  CEP  of  43.01  feet  slightly  higher  than 
the  CEP  of  41.22  feet  obtained  in  the  four  satellite  analysis, 
as  expected.  This  is  because  the  problem  is  essentially  planar 
and  three  satellites  are  sufficient  to  determine  position  and 
user  clock  bias  states.  The  barometric  altimeter  still 
stabilizes  the  INS  vertical  channel,  but  the  large  vertical 
position  errors  do  not  affect  the  planar  emitter  location 
problem.  Note  that  the  optimum  satellite  positions  in  the  true 
three  satellite  case  are  considerably  different  from  those 
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Table  4.7.  Three  Satellite  Final  Positions 


1 

1 

Azimuth  (-) 

Elevation  (-) 

LOS  Bias 

SV-1 

1 

1 

1 

299 . 6 

22. 4 

3 . 1 

SV-2 

1 

1 

1 

244 . 6 

11.8 

4 . 5 

SV-3 

1 

1 

105.9 

9 . 3 

3 . 7 

Table  4.8.  Three  Satellite  Error  Ellipse  Parameters 


1  Axis  1 

1 

(ft)  Axis  2 

Ancjle  from  Y  (-) 

CEP  1ft 

Obs-1 

1 

|  6.00 

I 

3 . 57 

62 . 8 

5 . 64 

Obs-2 

1 

|  5.99 

1 

3 . 61 

62 . 5 

5 .65 

Obs-3 

1 

|  5.51 

i 

3 .51 

59 . 8 

5.  31 

Emit 

1 

1  64.62 

8 .44 

-  0.12 

43.01 

obtained  in  the  pseudo-three  satellite  optimization.  The  true 
three  satellite  optimization  yields  satellite  positions  which 
correspond  to  those  anticipated  at  the  outset  of  this  study. 

SV-1  and  SV-3  have  assumed  positions  such  that  they  reduce 
observer  position  errors  primarily  in  the  east-west  direction 
and  have  also  essentially  aligned  themselves  so  that  LOS  bias 
errors  are  reduced.  This  can  be  seen  in  Table  4.7  where  SV-1 
and  SV-3  have  the  smaller  LOS  bias  terms.  SV-2  has  taken  up  a 
position  which  is  nearly  symmetric  with  SV-1  about  the  east-west 
axis . 


4.4  TWO  SATELLITE  OPTIMIZATION 


Again,  as  in  the  three  satellite  optimization,  only 
measurements  from  the  appropriate  number  of  satellites  are  used 
in  computing  the  error  covariance  matrix  updates.  This  presents 
a  problem  in  the  two  satellite  case  since,  in  general,  three  SVs 
are  required  to  solve  for  a  two-dimensional  position  and  user 
clock  bias.  For  the  two  satellite  case  in  this  study,  the 
navigation  filter  is  iterated  to  steady-state  using  measurements 
from  four  available  SVs  before  operation  is  degraded  to  two 
satellites.  This  allows  the  Kalman  filter  to  converge  on 
accurate  estimates  of  user  clock  bias  and  user  clock  drift 
states  before  degraded  operation  begins.  Using  these  state 
estimates  makes  it  possible  to  propagate  the  user  clock  bias 
forward  with  sufficient  accuracy  for  some  period  of  time. 


The  gradient  search  scheme  converges  on  the  optimum 
satellite  configuration  more  quickly  for  the  two  satellite  case 


than  in  the  searches  with  more  SVs.  This  is  due  in  part  to  the 
fact  that  both  satellites  move  onto  the  five  degree  constraint 
boundary  and  remain  there.  The  optimization  results  are 
presented  in  Fig.  4.11  and  Tables  4.9  and  4.10.  Satellite  paths 
to  the  optimum  positions  are  shown  in  Fig.  4.12. 

Note  that  for  the  two  satellite  case,  SV-2  has  again 
aligned  itself  along  the  east-west  axis  to  provide  maximum 
information  along  the  long  axis  of  the  emitter  error  ellipse. 
Additionally,  SV-1  has  assumed  a  position  that  provides 
essentially  the  same  amount  of  information  in  the  y  and  z 
directions.  The  increase  of  approximately  13  to  15  feet  in 
emitter  CEP  over  the  previous  cases  is  reasonable  since  the 
system  is  depending  upon  filter  estimates  for  user  clock  bias. 

It  is  significant  to  note,  however,  that  the  system  is  able  to 
achieve  good  performance  for  some  period  of  time,  even  when  only 
two  satellites  are  available  for  measurements. 

The  satellite  geometry  data  obtained  in  the  previous 
optimization  cases  is  now  used  to  develop  a  satellite  selection 
algorithm  that  will  minimize  emitter  CEP. 
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^.11.  Two  Satellite  Final  Geometry 
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Table  4.9.  Two  Satellite  Final  Positions 


1  Azimuth  (-) 

i 

Elevation  (-) 

LOS  Bias 

SV-1 

1 

|  143.9 

i 

5.0 

1.80 

SV-2 

I 

1  268.8 

5.0 

1.82 

Table  4.10.  Two  Satellite  Error  Ellipse  Parameters 


1 

Axis  1 

(ft)  Axis  2 

Anole  from  Y  (-) 

CEP  (ft) 

Obs-1 

1 

1 

1 

13.53 

2.96 

66.6 

9.70 

Obs-2 

1 

1 

1 

12.83 

2.98 

60.7 

9.30 

Obs-3 

I 

14.92 

2.97 

60.6 

10.54 

Emit 


82.95 


12.66 


0.89 


56.28 
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.5  SATELLITE  SELECTION  ALGORITHM 


When  developing  a  satellite  selection  algorithm,  there  is  a 
tradeoff  between  making  the  selection  criteria  restrictive 
enough  to  provide  acceptable  emitter  CEPs  and  liberal  enough  to 
allow  a  wide  range  of  possible  satellite  solution  geometries. 

The  basic  approach  used  in  this  study  is  to  examine  the 
sensitivity  of  the  emitter  CEP  to  variations  in  satellite 
position  from  the  optimum. 


The  previous  optimization  results  indicate  that  the  primary 
satellite (s)  should  be  aligned  along  the  semi-major  axis  of  the 
emitter  error  ellipse.  In  addition,  if  more  than  two  satellites 
are  available,  the  primary  satellites  tend  to  lie  along 
essentially  the  same  LOS.  Comparison  of  the  different 
optimization  results  is  useful  to  show  trends  in  solution 
geometry.  These  trends  are  very  useful  in  developing  selection 
criteria.  The  final  satellite  positions  for  the  four  satellite 
optimization  in  Fig.  4.5  show  SV-1  and  SV-2  very  close  to 
alignment  with  the  semi-major  axis  of  the  emitter  error  ellipse 
and  along  the  same  LOS  as  SV-3.  Note  that  SV-4  is  positioned  to 
provide  information  primarily  along  the  direction  of  the  semi¬ 
minor  axis  of  the  emitter  error  ellipse.  Comparing  these 
positions  with  the  pseudo-three  satellite  case  of  Fig.  4.7, 
where  SV-4  provides  no  horizontal  position  information, 
indicates  that  SV-1  and  SV-2  assumed  approximately  the  same 


positions,  while  SV-3  moved  to  a  position  that  provides  more 
information  along  the  semi-major  axis  to  compensate  for  the  loss 
of  information  from  SV-4.  In  fact,  SV-3  provides  equal 
information  in  each  axis  direction  of  the  emitter  error  ellipse. 

The  three  satellite  case  in  Fig.  4.9  indicates  SV-3  is  the 
primary  satellite,  while  SV-1  and  SV-2  have  taken  up  positions 
providing  information  primarily  along  the  semi-major  axis,  but 
also  some  information  along  the  semi-minor  axis.  It  is  the 
tradeoff  between  this  semi-major  and  semi-minor  axis  information 
that  is  instrumental  in  algorithm  development.  In  the  two 
satellite  case,  SV-2  is  the  primary  satellite  and  SV-1  moves  to 
the  equal  information  position  as  shown  in  Fig.  4.10.  These 
comparisons,  then,  form  the  basis  for  selection  criteria. 

Selection  criteria  are  summarized  as  follows: 

1.  Primary  satellite(s)  should  lie  along  the  semi-major 
axis  direction  of  the  emitter  error  ellipse  at  an  elevation 
of  10-25°. 

2.  Primary  satellites  should  lie  along  the  same  LOS  to 
reduce  LOS  bias  errors. 

3.  Additional  satellites  should  be  selected  to  provide 
adequate  information  in  the  general  semi-minor  axis 
direction. 


Satellites  are  selected  which  meet  the  criteria  or  come 
closest  to  satisfying  them.  Optimization  results  indicate  that 
as  the  primary  satellite (s)  move  closer  to  the  semi-major  axis 
direction,  the  additional  satellite (s)  move  toward  the  equal 
information  position.  The  three  satellite  case  illustrates  the 
situation  where  all  the  satellites  are  equally  dominant,  i.e. 
all  are  at  approximately  the  same  angle  with  respect  to  the 
semi-major  axis  of  the  emitter  error  ellipse.  The  selection 
algorithm  must  be  implemented  in  such  a  way  that,  as  less 
information  is  available  along  the  semi-major  axis  direction, 
additional  satellites  are  selected  to  provide  information 
primarily  along  that  direction. 

The  amount  of  change  in  emitter  CEP  as  the  satellites 
deviate  from  the  optimum  positions  is  addressed  somewhat  in  this 
study,  though  not  specifically.  Gradient  search  data  indicates 
that  large  changes  in  satellite  azimuth  can  be  made  with 
relatively  little  effect  on  the  CEP.  This  is  encouraging 
because  it  opens  a  wider  window  for  satellite  selection,  while 
maintaining  a  reasonable  CEP.  Data  regarding  CEP  change  versus 
satellite  elevation  is  not  so  readily  available,  but  changes  in 
elevation  in  the  range  of  interest  (up  to  15°)  should  not  result 
in  significant  increases  in  CEP. 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  CONCLUSIONS 

This  study  shows  that  emitter  location  errors  can  be 
significantly  reduced  by  selecting  satellites  according  to  the 
selection  criteria  presented  in  Chapter  4,  rather  than  by 
minimum  GDOP  criteria.  Three  satellite  performance  is  nearly  as 
good  as  that  obtained  using  four  satellites,  since  only  three 
satellites  are  necessary  to  determine  the  emitter  horizontal 
position  and  user  clock  bias.  When  operation  is  further 
degraded  to  the  two  satellite  case,  emitter  location  is  still 
slightly  better  than  that  obtained  using  four  satellites 
selected  to  minimize  GDOP. 

Although  the  global  minimum  of  the  CEP  cost  function  was 
not  found,  the  local  minima  obtained  for  all  of  the  fully 
determined  cases  agreed  very  closely.  This  indicates  that 
little  improvement  in  CEP  can  be  expected  beyond  what  is 
achieved  in  this  study. 

Iterating  the  error  covariance  matrix  to  steady-state  after 
each  satellite  movement  does  not  have  a  significant  effect  on 
either  the  final  optimized  satellite  positions  or  on  the  minimum 
CEP  obtained.  This  verifies  the  hypothesis  by  Lewantowicz  (4) 
that  using  only  a  single  iteration  of  the  navigation  Kalman 


filter  between  satellite  movements  without  reinitializing,  does 
not  significantly  affect  the  results.  Modifications  to  the 
gradient  search  routine  did,  however,  result  in  more  uniform, 
well-behaved  "tracks"  for  the  satellites  as  they  moved  to  their 
optimum  positions. 


5.2  RECOMMENDATIONS 


Since  the  weighted  gradient  search  routine  used  in  this 
study  is  highly  unlikely  to  converge  to  the  global  minimum, 
another  approach  should  be  used  that  has  a  higher  probability  of 
finding  the  global  minimum.  It  may  be  necessary  to  compute  the 
cost  at  a  very  large  number  of  random  points  in  the  eight¬ 
dimensional  Ps  space  to  identify  candidate  regions  in  which  to 
perform  gradient  searches  for  the  global  minimum. 

The  computational  loading  required  to  perform  the  weighted 
gradient  search  to  the  local  minimum  is  enormous.  This  is  due 
in  part  to  the  nature  of  the  gradient  search,  in  that  it  tends 
to  converge  at  a  rapid  rate  initially,  but  converges  very  slowly 
as  it  nears  the  minimum.  The  possilility  of  searching  initially 
using  the  weighted  gradient  method,  since  it  guarantees  that  at 
least  a  local  minimum  will  be  found,  and  then  switching  to  a 
more  rapidly  converging  algorithm,  such  as  the  Newton-Raphson 
method,  as  the  minimum  cost  is  approached,  should  be  explored. 
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The  selection  criteria  established  in  Chapter  4  should  be 
used  to  select  satellites  from  a  real-world  constellation  of 
available  satellites.  The  emitter  CEP  obtained  using  these 
"optimum"  satellites  could  then  be  compared  to  the  CEP  obtained 
using  satellites  selected  using  minimum  GDOP  criteria.  This 
would  provide  a  "real-world"  performance  evaluation  of  the 
emitter  location  system  proposed  in  this  study. 
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Appendix  A:  Cost  Gradient  Function  Derivation 

The  following  appendix  is  taken  in  its  entirety  from  the 
thesis  by  Lewantowicz  (4:104-107). 
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A.  1 


APPENDIX  3 


COST  GRADIENT  FUNCTION  DERIVATION 

The  cost  function  derived  in  Section  2.6  is 

J  -  tr  [EP+]  (B-1  ) 

where  E  is  an  n  *  n  constant  matrix  consisting  of  zeros  except  for  the 
two  diagor  1  elements,  corresponding  to  the  two  states  that  represent  the 
horizontal  emitter  position  errors,  which  are  =>  1.0.  The  matrix  P+  is 
the  n  *  n  symmetric  error  covariance  matrix  at  time  k.  However,  at  a 
fixed  k,  P+  is  a  function  of  only  the  variables  with  respect  to  which 
the  cost  J  will  be  minimized. 

The  positions  of  GPS  satellites  or  the  positions  of  observers,  or 
both,  could  be  the  variables  of  J.  In  fact,  Chapter  ^  covers 
optimization  with  respect  to  GPS  positions,  and  Chapter  5  covers 
optimization  with  respect  to  both  GPS  and  ooserver  positions.  Define  the 
vector  P3  as  the  vector  of  those  position  variables.  Then  the  cost  J 
is  a  scalar  function  of  P3. 

J  -  J  (P+( 

-  j  {p*[h]| 

-  J  (P+[H(Ps)]l  (B-2) 


where  H  is  an  m  *  n  measurement  sensitivity  matrix  and  P3  is  v.  i 
dimensioned  vector. 

In  finding  the  minimum  of  a  function,  tv-  -  — *  • 

in  several  minimum  cost  search  algorithms. 
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.  i  ---  . . . . . . . 


3. 


Using  the  chain  rule  of  partial  derivatives 


T  3J  3H 

*  3H  3P 

— s 


(B-4) 


where  the  first  term  is  an  m  *  n  matrix  and  the  second  term  is  m  *  n  *  2. 
third  order  tensor. 


The  first  partial  derivative  term  of  Equation  ( B-U )  is  derived. 


J  -  tr  [EP+] 


-  tr  [e(HTR"‘  H  +  P"  V*] 


Then  the  variation  in  J  is 


6 J  -  tr  [ESP+  ] 


where  6  is  the  variation  symbol 


(B-5) 


(B-6) 


6P+  -  -  (H 


TR“‘H  ♦  P"  ')  *  (  5HTR”  'H  ♦  hV1  6H)(HTR_1  H  ♦  P"  ) 


-  -  P*(6HT  R”‘  H  +  HT  R"'  6H)P  + 


(B-7) 


Therefore,  Equation  (B-6)  becomes 


however , 


iJ  -  tr  [-EP+(6HTR” 1  H  ♦  hV1  5H)P*] 
tr  [ab]  -  tr  [BA ] 

6J  -  tr[-P*EP*(6HTR*‘  H  ♦  HTR‘‘  6H )  ] 


(B-8) 


(8-9) 
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however,  tr[Aj  -  tr  [A  J 

Thus  the  variation  in  cost  J,  with  respect  to  the  measurement  sensitivity 
matrix  H,  becomes 

6J  -  tr[-P*EP*(HTR_,6H  *  HTR"'<5H)] 


-  -2  tr[p+EP+HTR“ l6H] 


(B-10) 


Using  the  result  attributed  to  Kleinman,  D.L.  in  [ 1 0 ] ;  given 


f(X)  -  tr[M(X>] 


<B-1 1 ) 


where  f  is  a  scalar  function  of  the  matrix  X,  and  M(X)  is  a  matrix  valued 
function  of  the  matrix  X  then  in 

f(X  ♦  eAX)  -  f(X)  -  etr[M(X)AX]  (B-12) 


as  e  >  0.  The  derivative  of  a  scalar  function  f  with  respect  to  the 
matrix  valued  variable  X  is 


3f(X) 

3X 


MT(X) 


(B-l 3) 


Thus  the  variation  Equation  (B-8),  expressed  as  a  partial  derivative 
becomes 


«  -  2[P*EP*H  R“ 1  ]  T 

on 

-  -  2R~'HP+EP+ 


{ B-1 4) 


where  3J/3H  is  an  m  *  n  matrix. 


Next,  the  3H/3£a  is  simply  an  element  by  element  partial 
derivative  of  H  with  respect  to  each  element  of  the  vector  P3.  For 
P3  of  dimension  l ,  there  are  l  m  *  n  matrices,  the  third-order  tensor, 
of  partial  derivatives.  Thus  each  component  of  the  cost  gradient  vector 
is  computed  according  to 


T  _  3J  3H 

£  3H  3P 

-s 


y  hit  i  j )  ( i  j ) 

J-1  -3 


( B-1 5 ) 


where  3J/3H  is  computed  as  shown  in  (B-H).  Computation  of  3H/3P3(k) 
matrices  is  presented  in  Appendices  C  and  D  for  GPS  only  and  for  GPS  with 
observer  optimizations  respectively. 
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